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Entanglement Perturbation Theory for Antiferromagnetic Heisenberg Spin Chains 



O 

(N 

D 

o 



u 

a 

o 
o 



00 

o\ 
o 



Lihua Wang^'^ and Sung Gong Chung^''^ 
^ Department of Physics and Nanotechnology Research and Computation Center, 
Western Michigan University, Kalamazoo, Michigan 49008, USA 
^ Computational Condensed-Matter Physics Laboratory, 
RIKEN ASI, Wako, Saitama 351-0198, Japan and 
^ Asia Pacific Center for Theoretical Physics, Pohang, Cyeonhuk 790-784-, South Korea 

(Dated: October 3, 2011) 

A recently developed numerical method, entanglement perturbation theory (EPT), is used to 
study the antiferromagnetic Heisenberg spin chains with z-axis anisotropy A and magnetic field 
B. To demonstrate the accuracy, we first apply EPT to the isotropic spin-i antiferromagnetic 
Heisenberg model, and find that EPT successfully reproduces the exact Bethe Ansatz results for 
the ground state energy, the local magnetization, and the spin correlation functions (Bethe ansatz 
result is available for the first 7 lattice separations). In particular, EPT confirms for the first time 
the asymptotic behavior of the spin correlation functions predicted by the conformal field theory, 
which realizes only for lattice separations larger than 1000. Next, turning on the z-axis anisotropy 
and the magnetic field, the 2-spin and 4-spin correlation functions are calculated, and the results are 
compared with those obtained by Bosonization and density matrix renormalization group methods. 
Finally, for the spin-1 antiferromagnetic Heisenberg model, the ground state phase diagram in A 
space is determined with help of the Roomany-Wyld RG finite-size-scaling. The results are in good 
agreement with those obtained by the level-spectroscopy method. 

PACS numbers: 71.10.Li . 75.10.Jm , 75.40.Mk 



I. INTRODUCTION 

The physics of strongly correlated systems has been 
one of the central interests in condensed matter physics 
since the discovery of high-temperature superconductiv- 
ity in 1986^. The successive discovery of a variety of 
exotic materials over the last two decades, showing su- 
perconductivity, magnetism and their co-existence, has 
further fuelled to our desire for the development of pow- 
erful numerical methods to study quantum many-body 
systems^. While the various kinds of mean field theory 
(MFT), the dynamical MFT and the density functional 
theory (DFT) being in the forefronlj^ii, have been ap- 
plied to many different systems to understand the strong 
correlation phenomena, there have also been continuing 
efforts concurrently for the development and refinement 
of non-MFT many-body methods. One of the promis- 
ing movements has its origin in density matrix renor- 
malization group (DMRG)^'^ and continues offering its 
refinements and modifications. Namely, an observation 
that DMRG leads to a matrix product representation of 
state vectors^, the infinite time evolving block decima- 
tion (iTEBD)*, and a variant using a variational prin- 
ciple rather than real space renormalization group (RG) 
to update local wave functions^, which advanced more 
recentljiiSiii to an improvement over the DMRG infinite 
algorithm^. The tensor extension in 2D began in the 
DMRG worki^, followed by many recent works, which 
primarily use RG to update local wave functions, includ- 
ing infinite projected entangled pair states (iPEPS) ^'^'^^ . 

The idea of yet another non-MFT, entanglement per- 
turbation theory (EPT) which we employ in this paper, 
was essentially put forward in the analysis of the clas- 
sical statistical mechanics, the 2D, 3D Ising models^^. 



There, the notion of infinite, translational invariant ma- 
trix, tensor product states, not via DMRG but from the 
viewpoint of singular value decomposition (see Section 
II), and particularly the matrix, tensor product repre- 
sentation of operators was introduced to fully implement 
the old, well-known tool of transfer matrix method. In 
other words, the translationally invariant matrix, tensor 
product representations of states and operators are the 
essential ingredient for EPT to be maximally efficient. A 
successive introduction of the notion of density matrix 
e~^^ with taking a mere parameter /3 — > 0, /3 is not the 
inverse temperature, reduces the quantum ground state 
problem to that of classical statistical mechanics^^J^. 

What are the differences among DMRG, its variants 
and EPT? They are all the same in that they are es- 
sentially variational, but they are all different in details, 
and details crucially matter here. All of these methods 
are after wave functions in the matrix, tensor product 
form, and the key step is how to update the local wave 
functions. Methods discussed in Refsi^i^iii^rii do it by 
real space RG. The variants^"^^ do it without real space 
RG but they keep the orthogonality condition on the lo- 
cal wave functions, a characteristic of DMRG, by a uni- 
tary transformation. In terms of our local wave func- 
tions, see Section II, it reads as = 1 or 
X)s=i '?*(^'')^ = 1- While the variants appear to improve 
on DMRG in ID, a real test would be in 2D where DMRG 
has not been quite successful, and where the orthogo- 
nality condition remains to be successfully implemented. 
While some progress has been made by^-^^^-, it is quite 
interesting to see how far one can go along these lines 
in 2D. On the other hand, EPT treats local wave func- 
tions with neither RG nor imposing the orthogonality 
condition both in ID and 2D. The penalty is a numerical 
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instability that sometimes occm'S in the iterative process. 
Now notice that EPT is based on a distinctive philosophy, 
namely, the idea of divide and conquer realized mathe- 
matically by singular value decomposition (SVD). SVD 
allows us to write everything from state vectors to op- 
erators in the matrix, tensor products. Besides, there is 
no reason for the orthogonality condition in EPT. There- 
fore, EPT is undoubtedly simplest, and thus has a great 
potential for 2D only if we can establish a numerically 
stable algorithm. So far, we have succeeded in ID, as 
will be demonstrated in this paper, quasi- ID in another 
paper-l^, but only partially yet in 2D^^'^^. 

The model system we study in this paper is the ID 
antiferromagnetic (AF) Heisenberg model described by 
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, where = S*^ ± iS'^ are the spin flip operators. J is 
the exchange coupling (taken to be 1) and A denotes the 
exchange anisotropy. A = corresponds to the xy model, 
1 the XXX model, oo the Ising model and general A the 
xxz model. Despite its simplicity, the Heisenberg spin 
chain fascinated physicists for generations. It is now one 
of few many-body quantum systems that are understood 
fully, by both analytic and numerical tools. 

The first analytic breakthrough is Bethe ansatz 
(BA)^"— . It rigorously solves for the ground state (GS) 
and the excitation spectrum, though a complete equa- 
tion is computationally challenging to solve^i^— . Some 
rigorous results of BA are essential ingredients for other 
theories developed later such as bosonization^I and con- 
formal field theory (CFT)^. Nevertheless BA has lim- 
itation to calculate the long-distance correlations^. In 
contrast, those effective field theories emphasize on the 
long-distance behavior of the model. 

Bosonozation transforms the xxz model into a quan- 
tum sine-Gordon model expressed in bosonic field 
operatorsSfl. In the xxx model, the so-called Umk- 
lapp scattering term , from S^S^ interaction, vanishes 
marginally after iterative renormalization group trans- 
formations, with its fixed point being the Gaussian 
Hamiltonian2£i2i. Due to a spin rotational symmetry, it 
is then shown^ir— that the spin-spin correlation decays 
as , r being the spin separation. 

The GET, on the other hand, shows thatj^ there is a 
multiplicative correction, (Inr) = , to the spin-spin corre- 
lation functions, through a renormalization group analy- 
sis of the conformal dimension which receives correction 
from the (marginally irrelevant) Umklapp term. The ex- 
act long-distance behavior {2Tr)~^ {Inr)^ is derived 
in"^^. The power of GET lies in the highly stringent 
conformal invariance in the theory of complex functions, 
leading to the determination of the entire spectrum solely 
in terms of the central charge and confromal dimensions, 
the latter in turn leading to the spin-spin correlation 



functions in the long distance^. 

On the numerical side, there are a variety of works. 
The typical ones are Quantum Monte Garlo^i^ and 
DMRG^i22ii2S. Recently researchers are numerically ex- 
ploring spin chains with impurities'**', the mixed spin-1- 
spin-i chain^* and the spin chain with anisotropic cou- 
pling strength in the x, y, z directions^. In another 
EPT work, the EPT-e algorithm^?, following Eeynman's 
line^, was reported for the excitation spectrum of the 
spin chains, in particular the spin triplet excitation spec- 
trum of the spin-1 chain for the whole Brillouin zone 
including the Haldane gap aX k — ii. EPT-e calculations 
rely on EPT-g calculations to get a very accurate ground 
state wave function. 

After describing the detailed formalism in section II, 
we apply EPT to study the spin i AE Heisenberg chains 
in section HI, where we first demonstrate the accuracy of 
EPT in comparison with BA results in Sec.IIIA, and then 
focus on a macroscopically long-distance behavior of the 
spin-spin correlations, and thus examine how far one can 
go with the size of entanglement p, roughly a measure of 
correlations. In this regard, we could go as far asp = 209 
(main time consumption in this calculation is taken by 
an eigenvalue decomposition of matrices with rank of 
and a generalized eigenvalue equation of matrices with 
rank of 2p^ . See the definition of the rank of matri- 
ces in section II) for the spin-i xxx model. In Sec.IIIB, 
Bosonization's prediction on the anisotropic Heisenberg 
chain in a finite external magnetic field is confirmed. In 
particular, EPT has turned out to be more accurate than 
DMRG in the 4-spin correlation functions. However the 
critical exponent given by Bosonization for the xxx model 
is found to be slightly different from the exact result. 
Thus another field thoery, GET, is examined in Sec.IIIG. 
Its prediction on the asymptotic behavior, which is found 
to start around 1000 in lattice separations, is confirmed 
for the first time to our best knowledge. Next, in sec- 
tion IV, we study the ground state phase diagram of the 
spin-1 xxz AE Heisenberg chain, combining EPT with 
the Roomany-Wyld RG-finite-size-scaling. We can still 
use EPT-g algorithm, that is originally designed for the 
ground state, to look for the first excited state of this 
spin chain. This is based on the fact that the first excited 
state is at the wave-number k = tt, which can be han- 
dled by the bipartite local wave functions ^^'^ employed 
in EPT-g. Our result agrees with the phase diagram ob- 
tained by the level-spectroscopy method^^. Finally, we 
conclude the paper with discussions. 



II. EPT 

EPT first introduces the local wave functions and then 
couples them by the entanglement to represent the trans- 
lationally symmetric wave function. Second it reaches a 
general yet simple algebraic procedure to solve for the 
local wave functions, during which the uniform EPT will 
utilize the translational symmetry of the Hamiltonian to 
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FIG. 1: Schematic figure of the wave function 



greatly simplify the formulation. The first step is com- 
mon to all the EPT algorithms, therefore we will dis- 
cuss it first. In fact, the resultant system wave function 
in EPT expressed by a successive product of the local 
wave functions precisely falls into the matrix product 
state concept^. As for the second step, we have differ- 
ent ways to handle the operator, i.e., we either deal with 
the Hamiltonian itself or the density matrix, e~^^ . It is 
how to efficiently find the local wave functions that dis- 
tinguishes EPT from other MPS methods^^^^i^^iS. Be- 
cause at the beginning of EPT development it was the 
second method, i.e. dealing with the density matrix, that 
was applied to the 2D&3D Ising models, it is called EPT- 
gl. The EPT-g2 algorithm was developed later to handle 
the Hamiltonian directly. In fact, EPT-g2 formulation is 
more lengthy. So we will discuss EPT-gl first in the fol- 
lowing SecinS and then EPT-g2 in SecHTBl The com- 
parison between them is given in Sec lII CI 

We now take a 4-sites spin chain as an example to il- 
lustrate how to write the system wave function in terms 
of local wave functions. The number of local spin states 
is denoted as g. So 5 is 2 for spin-i and 3 for spin- 
1. The wave function of this spin chain is written 
as V' (oi, 02, 03, 04), which is a q'^ dimensional vector. 
However it can also be regarded as a q x matrix 
"0(01,020304). We then use singular value decomposi- 
tion (SVD) to separate oi from the rest as 



V' (01,020304) = Aa^^a ^aBa2 

^ A' R 

^ai,a-^a2a3a4,a 

where and below the repeating index means summation. 
The eigenvalue Aq is absorbed into Aa^^a to give A'^^ ^. 



We then SVD -6020304,0 to separate 02 from 0304, we have 

-^020304,0 — C a2 ^ah^lhD a^a^^b 

We keep using SVD until all the lattice sites are sepa- 
rated from each other and the eigenvalues are properly 
absorbed into the corresponding left matrices, which can 
now be regarded as the wave functions of local states 
01,02,03,04. Finally the system wave function ends up 
as 

ip (01,02,03,04) Ca2,ab^a3, 60-^04, c (2) 

The above shows how to use SVD to rewrite the wave 
function as a successive product of local ones. We bor- 
row the terminology 'entanglement' from quantum infor- 
mation theory to denote the coupling indices a, b, c, d in 
equation Now recall that the spin chain is transla- 
tional symmetric and has a periodic boundary condition. 
We should introduce some uniform local wave function 
^ to represent the vectors A',C',E',F. Considering the 
possible or |t spin arrangement of AF spin chains, we 
need two local wave functions and f ^ to include such a 
bipartite structure. Therefore we have the SVD-ed wave 
function for a general spin chain as follows 

^ = • • • Cba, » CL2 ® ^cdas ® Cdea, " ' " (3) 

where, the subscript 1 on the shoulder refers to the odd 
sites of the chain while 2 the even sites. The indices 
{01,02,03,04} account for the local spin states, running 
from 1 to 2 for spin-i, 1 to 3 for spin-1. The indices 
{a, b, c, d, e} are the entanglement indices, which are in- 
tegers, uniformly coupling the local wave functions over 
all spin sites. Fig H] illustrates how the system wave func- 
tion is expressed as a successive product of the local wave 
functions. 

Given the system wave function in the form of ([3l). 
we show how to derive the main equation to solve for 

and in general. First recall the way to derive the 
Schrodinger equation by the variational principle, where 
one varies the expectation value of the Hamiltonian (en- 
ergy) with respect to the system wave function. In the 
algebraic form, the Schrodinger equation is an eigenvalue 
problem. One now has to similarly vary the expectation 
value of a functional of the Hamiltonian with respect to, 
for example, as follows. 

± ( {^ie,e)\nH)\i^ie,e)) \ ^ „ 

We thus arrive at a generalized eigenvalue equation, 

x^{e,e)e = f^Y4e,e)e * = i,2 (s) 

The GS of the system is corresponding to either the 
smallest or the largest eigenvalue of ([S]), depending on 
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the functional f{H). We see that it is possible to make 
([5]) self-consistent by starting with two vectors and 
as a seed, and to determine Xi and X2, hence to deter- 
mine a set of S,l and as the new seed. The iteration 
is continued this way until the convergence is reached. 
Note that the size of the generalized eigenvalue problem 
is controlled by the entanglement, p. More precisely, it 
is 2p^ for spin-i and for spin-1. The larger entan- 
glement will give the more precise representation of the 
wave function. If it is large enough, the wave function 
could be exact. So we start with a small entanglement, 
say 1, to get the converged and for it. Then we in- 
crease the entanglement and solve for the corresponding 
^'s. We compare the physical quantities calculated from 
the local wave functions with increasing entanglements 
to seek for the convergence. EPT calculations show that 
the entanglement necessary for the convergence increases 
slower than a linear function of the system size and that 
the convergence for each entanglement is rather quick. 
Therefore it has high efficiency, which enables one to cal- 
culate very large system precisely because now the size 
to be handled is roughly of order vs. in the ex- 
act diagonalization. In short, EPT is a method of divide 
and conquer realized by SVD, Eq.(3) — divide and the 
generalized eigenvalue equation, Eq.(5) — conquer. 



A. EPT-gl 



Since any eigenstate | ipi) of H is also the eigenstate 
of density matrix e~^^ , the GS now bears the largest 
eigenvalue e~^^^ . The eigenvalue problem we try to solve 
in EPT-gl is, 

e'f"" I = e-^^' I (6) 

Let us consider for an illustration purpose the xxx 
Heisenberg Hamiltonian, and express it as the sum of the 
local bond Hamiltonian between two nearest neighbors 



H — ^ Hbond 
bond 



(7) 



with 



Hbond = {S^S^_^_-^ + S^^^) + JS-S-^-^ (8) 

We then choose an appropriate small positive /3 (0 < /3 ^ 
1) to safely separate the even and odd sub lattices in the 
density matrix 



e - = lie 

bond 



Now let us linearize e ^^>"'^^<' , again, due to the fact that 

0<;3< 1 



-I3H,, 



(10) 



where ^la takes four operators 1, , Sj^ , Sf on site i and 
Qa likewise operators on site i -I- 1. The index a runs 
from 1 to 4. The local density matrix can be written as 



{Ik I e-l^"^"--' I ij) w /„,,fe (8) 
where, for a spin-i chain. 




(11) 

(12) 
(13) 
(14) 



and so on. For a spin-1 chain, /'s and 5's should be 3 x 
3 matrices. 

Now let us write down the matrix representation of the 
even and odd bonds in the density matrix as follows 

• • • fa®ga ® f-^®g-^■■■ 
■■■fl3®9|3® fs^gs--- 

The vertical alignment above means that the matrix rep- 
resentations are at the same site. The whole density ma- 
trix will be written as 



K 



= •r 



-Hi'" 



(15) 



Fig 12] illustrates the expectation value of the density ma- 
trix. 

K is now substituted into Q as the functional / to 
have 

{i^{e,e)\K\i^{e.e))=Tr(A^i^) (i6) 
{^{e,e)\^{e,e))=Tr(B^/^) (i?) 

where the matrices A and B are the translationally sym- 
metryic units illustrated in Fig|3l Explicitly 



A-mi,ni — £,aeai^eca2^bfbi^fdb2^aP,ata2^Pt,bib2 
Bm,n = ^aeai^eca2^bfai^fda2 

mi = (7 — 1) X + (6 — 1) X p + a 
rii = (a — 1) X p'^ + (d — 1) X p + c 
m~{b — l)xp + a 
n = {d — 1) X p + c 



(18) 
(19) 
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FIG. 2: Schematic figure of expectation value of the density 
matrix in EPT-gl 
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trace. We then diagonalize the matrices A and B: 

A =RAL'f (21) 
B ^R'A'L"' (22) 

where R, A, and L are respectively the right eigenvector 
matrix, the diagonal eigenvalue matrix and the left eigen- 
vector matrix. Equation pop is then further rewritten as 
follows 



(Ai)^/"-^ (M) i?, =a^(A;//' 'mSB)R'^ (23) 

where i runs from 1 to Ap^ and j runs from 1 to . Now 
we can explicitly work out the variation with respect to 
the local wave functions, for example, Using equa- 
tion (|18l) . ((23)) leads to the detailed form of one of the 
generalized eigenvalue equation ([5]), 

Xi{aeai,bfa2)^lfa^_ = ^iYi{aeai,hfa2)S,lfa^ (24) 

where 

Xi(aeai, 6/02) 

_ (^2 f2 pi p2 p r (K \L/2-l 



Yi{aeai,bfa2) 



(25) 



(26) 



Likewise, we get the generalized eigenvalue equation for 
We solve for and as stated in Sec|lTl We have 
next step ^^■^ and repeat calculation until convergence, 

^1,2 „ fl^2 




FIG. 3: Schematic figure of forming the translationally sym- 
metric units in EPT-gl 



Above a,b,c,d,e, f run from 1 to p, 01,02,61,62 I'un 
from 1 to 5 and a, 13, 7 run from 1 to 4. We repeat that 
p is the entanglement and (7 is 2 for spin-^, 3 for spin-1. 
We see that ^ is a Ap^ x Ap^ matrix for spin-i and i? is a 
p^ X p^ matrix for both spin-^ and spin-1. The variation 
(HI) leads to 



(20) 



B. EPT-g2 

EPT-g2 handles the Hamiltonian directly. The GS en- 
ergy and first excited state energy can be solved by look- 
ing for the first and second minimum of the quantity 



(V' I g I V') 



(27) 



Now since both the wave function and the Hamiltonian 
are translationally symmetric, we rewrite (l?7)) as follows 



_ L I H§^:i, I ^) L (V. I H^Zl I ^) 
2 I ^) "^2 (V- I ^) 

We have the matrix representation 

^bib2Mia2 — (61 I (62 I Hbond \ Al) | 0,2) 



(28) 



(29) 



The symbol * denotes matrix product and then take a 



where | ai) and | 02) refer to the local basis vectors on two 
neighboring sites while (61 | and (62 | are the conjugate 
counter parts. Now we are allowed to write ([^5)) further 
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FIG. 4: Schematic figure of the formation of matrixes A, B, 
C and D in EPT-g2 



Above a,6, c, d, e, / run from 1 to p and ai, 02, 61,62 
run from 1 to q. We can write down the matrix C and 
D likewise. FigU] schematically shows the formation of 
matrixes A, B, C and D. If we vary ([50)1 with respect to 

and we will arrive at 

L/2-1 

L/2-1 

(33) 



We then diagonalize the matrices A and C: 

A =RAL'' 
C =R'K'L'^ 
Equation (155)) is rewritten as follows 



(34) 
(35) 



as follows 



LTr{A^/^-^B) LTr{C^/^-^D) 



2 Tr(A^/2 



2 rr(A^/2 



where 



A — tl t2 ^1 f2 
^m,n — SaeaiSeca2''6/aiS>/da2 

Bm.n — S,aeai^eca2^bfbi^fdb2^bib2.aic 

m=(b — l)xp + a 
n — (d — l)xp + c 



(30) 



(31) 
(32) 



e^(A,)^/'"'L. (M) Ri 



i=l 

p2 



i=l 

h{L,R,K,5A)+J^{L',R',K',5C) 



(36) 



where 



A(L,i?,A,M) 
f^{L',R',k',6C) 



p2 

E 

p2 

E 

m,n— l;m^n 

LmBR, 



A _ a'-^ 

i I777, i Lyj 



A' -A' 

m n 



{SA)R 

m 

m— 1 

p2 

i^,„„L; (,5C) + (L/2 - 1) ^ A;^/2~2i.,_i;^ (jc) R'„ 



Like in EPT-gl, we solve ([S]) and hence the problem. 
Note that the smallest eigenvalue and corresponding 
eigenvector give the GS energy and the local wave func- 
tions. The second smallest will bring about the first ex- 
cited state energy and the corresponding wave functions 
as long as the excited state has an appropriate wave func- 
tion structure. 



C. Summary 



In the last two subsections, we saw that both EPT- 
gl and EPT-g2 can solve for the GS. The EPT-gl is 
simpler, but the calculation is heavier because the ma- 
trix A for EPT-gl is larger than B and D for EPT-g2. 
However, EPT-gl is especially suitable for infinite sys- 



7 



chain length 


EPT rcsuh 


Bcthc Ansatz result 


16 


-0.4463935 


-0.4463935 


64 


-0.4433459 


-0.4433485 


256 


-0.4431551 


-0.4431597 



TABLE I: Comparison of the GS energies between Bethe 
Ansatz and EPT 



terns with translation symmetry, because only the largest 
eigenvalue of A needs be retained. On the other hand, 
EPT-g2 has matrices E and F containing all the eigen- 
values and eigenvectors of A, B, C, D, hence it becomes 
slower than EPT-gl for infinite systems. Actually we 
used both EPT-gl and EPT-g2 (/3 is set to be 10"*^ for 
up to thousand spin sites) for finite spin chains. Since 
EPT-g2 does not make any approximation and the two 
methods give exactly the same results and converge with 
the entanglement at the same speed, the small parameter 
/? in EPT-gl algorithm is not an issue. Later we will see 
EPT-gl calculation for infinite spin-i chains gives very 
accurate results. Moreover EPT-g's can be applied to 
the first excited state for spin-1 chains as well, because 
the first excitation occurs when the wave number k — tt 
which means the bipartite trial wave function is suitable. 
On the other hand, the first excited state in spin-^ chains 
cannot be calculated by EPT-g's since it does not have 
this character. We then need to use EPT-e as discussed 



III. SPIN-i CHAINS 

First we use EPT-g to solve for the GS without the 
magnetic field, and calculate the energy, the local magne- 
tization and spin-spin correlations. We especially investi- 
gate the long-range behavior of 2-spin correlations, con- 
sidering a very large entanglement and the infinite size, 
and examine the prediction of CFT^^. As mentioned in 
Sec HTCl EPT-g2 works as well as EPT-gl for finite sys- 
tems. But for infinite systems, EPT-gl is more suitable. 
Second we calculate the model with z-axis anisotropy and 
in the external magnetic field. In particular, we stud- 
ied both the long-range 2-spin and 4-spin correlations, to 
check the predictions of the bosonization method^. 

One of the merits of EPT is that it can give precise 
local wave functions from which many quantities can be 
calculated. However, one must ask how to judge the local 
wave function's convergence. We can choose one of those 
quantities as a convergence indicator. The convergence 
speeds of the GS energy and the local magnetization for a 
256-sites chain are shown respectively in Figl5]and Fig|6l 
We see that the convergence of the GS energy is reached 
before entanglement p = 25 and the local magnetization 
is converged after p = 30. FiglT] shows the convergence 
of another important quantity for the same chain: the 
spin-spin correlations, Wa{r) = {SqS'^). It is shown that 
Wz and Wx curves start to coincide after p = 35, which 
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FIG. 5: Convergence with entanglement of the GS energy for 
a 256-sites spin-i chain 
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FIG. 6: Convergence with entanglement of the local magne- 
tization for a 256-sites spin-i chain 



is a sign of the convergence in 2-spin correlations because 
the symmetry among x, y and z axis are recovered. In 
short, the convergence of the correlation functions ap- 
pears to be the ultimate indicator for the convergence of 
all quantities. 



A. Comparison with Bethe ansatz 

We have calculated the spin-i chains of length 16, 
32, 64, 128, 256 and 512. EPT agrees well with Bethe 
ansatz^- for the GS energies, as shown in TabU] 
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FIG. 7: Convergence with entanglement of the spin-spin cor- 
relations for a 256-sites spin-i chain 



FIG. 8: Log-log plot of the spin-spin correlation functions over 
the first 501 sites for an infinite spin-i chain at entanglement 
180 



distance 


EPT 256-site 


EPT infinite 


Bethe Ansatz 


1 


-0.1477158 


-0.1476404 


-0.1477157 


2 


0.0606715 


0.0606355 


0.0606798 


3 


-0.0502196 


-0.0501772 


-0.0502486 


4 


0.0346027 


0.0346248 


0.0346528 


5 


-0.0308457 


-0.0308430 


-0.0308904 


6 


0.0243621 


0.0244263 


0.0244467 


7 


-0.0223980 


-0.0224636 


-0.0224982 



TABLE II: Comparison of the spin-spin correlations between 
EPT and Bethe Ansatz over the first 7 sites 



To check the accuracy of the calculated correlation 
functions, we compare in TabHIl EPT with the Bethe 
ansatz (generating function method)^^ for the first 7 lat- 
tice separations, which are the only spin-spin correlations 
handled by Bethe ansatz up to date. We see that EPT 
agrees with BA. Note that at entanglement 180, in the 
third column of Tab|lTl EPT for an infinite chain achieves 
the relative error of 0.1% throughout the first 7 lattice 
separations. 



B. Comparison with Bosonization 



better prediction with a multiplicative logarithmic cor- 
rection was given by CFT, which will be discussed in 
Sec.IIIC. 

For the moment, let us consider the Bosonization's pre- 
diction on the xxz model in external magnetic field B, 
where the bosonization makes precise prediction for both 
2-point and 4-point functions since the asymptotic loga- 
rithmic correction is absent for non-zero B. The 2-spin 
correlation function is predicted to be^^, 



(5f 5f+,,) =M' 



1 



n=0 



i(-ir 



cos [{2n + 1) Qr] 

^(2n+l)V') 



(37) 



where the decay exponent is a function of A and the 

magnetization Af , — -f- and Q ~ 27rAf is an incom- 
mensurate wave number. The 4-spin correlation func- 
tions are likewise given by^^, 

(■ i^^^l+l + ^l+l) i^l+r^l+r + l + ^l+r^l+r+l) ') 



(-1)'' 

.16Bi^^-^cos(Qr) 



165^ 



(38) 



Let us compare Bosonization's prediction on the spin- 
spin correlations with our result. In fig El the linear 
fit of the log-log plot of the spin-spin correlations for 
odd separations shows a power-law decay, Wz{r) ~ 
-0.1473r-0-9604 ^ith a fitting error of 0.1%. This fit- 
ting also applies to the even r but with a positive sign 
reflecting the AF nature. We see that the critical ex- 
ponent given by EPT is —0.9604, indeed very close to 
Bosonization's prediction of —1 but slightly different. A 



(39) 



where : 3? := 3? — (5R). In principle, the above coeffi- 
cients Bi, B'q, B2, a„ and so on can be determined by fit- 
ting to numerically computed correlation functions. Au- 
thors o:^ indeed used DMRG to determine these quan- 
tities. We have calculated the two and four spin corre- 
lation functions for the anisotropy A — 0.5, 0, —0.5, and 
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B 


A = -0.5 


A = 


A = 0.5 


M=0.05 


0.03395 


0.07820 


0.13480 


M=0.25 


0.16058 


0.35350 


0.56680 



TABLE III: The magnetic field intensities for difi'erent com- 
binations of M and A. 
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FIG. 9: M converges to 0.05 suddenly at entanglement 24. 
B = 0.1348 and A = 0.5 



the magnetization M = 0.05,0.25 for a 200-sites spin-^ 
chain, as'^^ did. The difference is that we used the pe- 
riodic boundary condition while^ used the open bound- 
ary condition which is more accurate than the periodic 
boundary condition in the DMRG calculation. We vary 
the magnetic field B to gradually approach an aimed M 
for each given A. Tab lllll gives the magnetic field inten- 
sities for some combinations of M and A. In fact, one 
can vary B to have the calculated M close to the aimed 
value. This can be done quickly at small entanglement. 
Then a convergence is checked with respect to increasing 
entanglement. One sees that the calculated M suddenly 
converges to the aimed value in FiglQl For the 2-spin 
correlation, we mainly compare (SqS^) with (b) of Fig. 1. 
irt^, while {SqS!^) is quite similar to (a) of that figure 
(we use the periodic boundary condition). Fig llOl shows 
that {S§S^} by EFT, unlike DMRG, does not have the 
artificial oscillation at the edges. 

We now present in Fig llll the 4-spin correlation func- 
tions by EFT for A = -0.5,0,0.5 and M = 0.05. 
The calculations are all converged with the entangle- 
ment. A small deviation near the center of the chain 
from the power-law decay is due to the periodic bound- 
ary condition. {S^ Si^-j^Si^_^.^Si^_^.^^-j) and the envelope of 
(: {S^Si_^_^+ Si Si[^_^) :: (5*;+ ^S*;,^^,^^ -f 5;^^^;+ J :} 
decay in a power-law as a linear fitting in a log-log plot 
shows, which is predicted by'^^. However, we see a clear 
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FIG. 10: {SgSr) by EPT does not iiave the artificial osciUa- 
tion at the edges. 



period-10 oscillation in plots (a) and (6) while DMRG 
does not show this (Fig. 3. ini^). In fact, plot (c) also has 
an oscillation with the period of 10, although subtract- 
ing the mean value from the correlation makes it hardly 
visible. This oscillation was clearly seen for the first time 
owing to EPT's great precision. Note that Q — O.Itt 
when M = 0.05 and the absolute values were taken in 
the plots. It is exactly what was predicted by the term 
cos{Qr) in 



C. Comparison with CFT 

Because of a required high precision, CFT's prediction 
of logrithmic correction to a power law decay of spin- 
spin correlations has not been addressed by any methods 
so far, to our best knowledge. In Fig. 12, we plot Wz 
and confirmed the CFT result. Note that the conver- 
gence of Wz with entanglement to CFT is seen around 
r ^ 1000 for the entanglement p > 180. This also shows 
that an asymptotic behavior is realized for r 1000. 
At this extreme accuracy, the obtained ground state en- 
ergy -0.4431467 is very close to the exact value of Bethe 
ansatz -0.4431471. It would be interesting to check an- 
other prediction of CFT^^ and series expansion^": if and 
how the asymptotic correlation functions exhibits critical 
behavior at A = 1 



IV. SPIN-1 CHAINS 

Similarly we calculated the spin-1 chains of many 
different lengths. We compare the EPT results with 
DMRG. The GS energy for a 48-sites chain is -1.401482 
by EPT at entanglement p — 20 versus -1.401484 by 
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(a) A = 0.5, Af = 0.05, L = 200 




(b) A = 0,M = 0.05, L = 200 
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(c) A = -0.5, M = 0.05, L = 200 
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FIG. 11: The four-spin correlation functions for M 
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FIG. 12; Spin-spin correlation functions for an infinite spin-i 
chain 



0.0 
-1.5 ■ 
-3.0 
g -4.5 ■ 
-6.0 



o 



-7.5 



-9.0 ■ 




-i — I — I — I — I — I — I — I — I — I — I — I — I — I — I 
10 20 30 40 50 60 70 



FIG. 13: The spin-spin correlation functions exponentially 
decay with the distance 



DMRG^. The energy gap between the GS and the first 
excited state in the same calculation is 0.41242 by EPT 
and 0.41232 by DMRG. The spin-spin correlations were 
calculated from the GS wave function. We see a clear 
convergence in the semi-log curve of the correlation func- 
tions even as early as at p = 18 for a 64-sites chain in 
FiglTSl The straight line in the semi-log plot indicates 
an exponential decay of the spin-spin correlation with 
the distance, consistent with the existence of the energy 
gap in spin-1 chains. 

We have also calculated the GS phase diagram of the 
spin-1 xxz chain, A ^ 1, using the Roomany-Wyld RG 
finite-size scaling-^. In FiglHl the four phases are iden- 



11 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 



I 

FIG. 14: Phase diagram of the spin-1 xxz chain 



tified in the plot of Gap(L) x L, where Gap(L) denotes 
the energy gap as a function of the chain size L. When A 
is fairly large, it approaches the Ising limit. The GS is a 
unique Neel state. The spin-1 xxz chain enters from Neel 
phase to Haldane phase at A = 1.192, where the energy 
gap disappears. From this point to A = there is an en- 
ergy gap for all chain length. And at A = 1, the energy 
gap converged to a finite value, the Haldane gap. From 
A = OtoA = — lis the XY phase, the infinite chain has 
no energy gap. At A = —1, all finite chains have no en- 
ergy gap. After this point, it enters ferromagnetic phase, 
where the GS is doubly degenerate. Our phase diagram 
agrees with that of the level-spectroscopy method^. 

V. CONCLUSION 

Two EPT algorithms, EPT-gl and EPT-g2, have been 
applied to calculate the ground sate properties of AF 



spin chains. They give exactly the same results with 
the same convergence speed. Among them, EPT-gl is 
especially suitable for the infinite system. By successfully 
comparing with Bethe ansatz, CFT and Bosonization, we 
have seen the EPT's ability to handle quantum systems 
of arbitrary size. This is especially a good news to field 
theory because some of their predictions have not been 
verified by any numerical calculations before. Along with 
its further ability to calculate excited states either by 
EPT-g algorithm or EPT-e^, EPT's broad apphcability 
in ID quantum systems will be clear. 

Among many possible applications and desirable ex- 
tensions of EPT, we are currently examining the non- 
uniform matrix product states in the spin-^ AF Heisen- 
berg models in ID and 2D. In particular, the Ji — J2 
model in 2D including the square and triangular lat- 
tices. Since such a non-uniform system is more or less 
inherently of short-range nature, one would not expect 
much improvement over DMRG. A question is rather if 
the non-uniform EPT can still maintain its simplicity 
and stability in algorithm and numerical efficiency. A 
work is currently underway, and will be reported in the 
near future. A true new aspect, however, might emerge 
when we consider impurities/nano-structures embedded 
in strongly correlated host materials of macroscopic ex- 
tension, possibly realized in Kondo phenomena. It is 
quite interesting to see how EPT works here. 
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